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ABSTRACT 

We show how observations of multiply-imaged quasars at high redshift can be used as a probe 
of dark matter clumps (subhalos with masses 10 9 M©) within the virialized extent of more 
massive lensing halos. A large abundance of such satellites is predicted by numerical simulations 
of galaxy formation in cold dark matter (CDM) cosmogonies. Small-scale structure within galaxy 
halos affects the flux ratios of the images without appreciably changing their positions. We use 
numerical simulations to quantify the effect of dark matter substructure on the distribution 
of magnifications, and find that the magnification ratio of a typical image pair will deviate 
significantly from the value predicted by a smooth lensing potential if, near the Einstein radius, 
only a few percent of the lens surface density is contained in subhalos. The angular size of the 
continuum source dictates the range of subclump masses that can have a detectable effect: to 
avoid confusion with gravitational microlensing caused by stars in the lens galaxy, the background 
source must be larger than the optical continuum-emitting region of a QSO. We also find that 
substructure will cause distortions to images on milli-arcsecond scales and bias the distribution 
of QSO magnification ratios - two other possible methods of detection. 

Subject headings: cosmology: theory - dark matter - galaxies: formation - gravitational lensing 

1. Introduction 

The popular model of hierarchical structure formation in a universe dominated by cold dark matter 
(CDM), while quite successful in matching the observed large-scale density distribution, is currently facing 
a 'small-scale crisis' (e.g., Moore 2001). First and foremost, CDM simulations appear to produce halos 
that are too centrally concentrated compared to the mass distribution inferred from the rotation curves of 
(dark matter-dominated) dwarf galaxies (Moore et al. 1999b; Klypin et al. 2001; Navarro, Frenk, & White 
1997). In addition to the cuspy density profiles, the predicted abundance of satellites of mass greater than 
10 8 M Q within the virialized extent of a Milky Way's halo is more than an order of magnitude larger than 
the number of dwarf galaxies with comparable mass observed within the Local Group (Moore et al. 1999a; 
Klypin et al. 1999b; Mateo 1998). A number of related but more complex problems is also being discussed 
in the literature (see Sellwood & Kosowsky 2001 for a recent review), together with the significance of some 
of the observations (van den Bosch & Swaters 2001). It is unclear at this stage whether the CDM paradigm 
is actually incorrect and needs to be modified to include, e.g., self-interacting (Spergel & Steinhardt 2000) or 
warm dark matter (Bode, Ostriker, & Turok 2000; Colin, Avila-Reese, & Valenzuela 2000), or whether some 
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of these discrepancies may have a more 'astrophysical' origin. Feedback processes such as photoionization, 
for example, may prevent gas cooling and inhibit star formation in the majority of the dwarf subhalos (e.g. 
Bullock, Kravtsov, & Weinberg 2001). In this case the dark matter satellites of the Milky Way would still 
be present, but only a small fraction of them would actually be visible, having formed stars prior to the 
reheating of the intergalactic medium. 

Gravitational lensing is a powerful tool for studying the structure of dark matter halos. The statistics 
of wide-separation lenses can be used to constrain the amount of mass in the cores of dark matter halos on 
group and cluster mass scales (Keeton & Madau 2001; Li & Ostriker 2001; Flores & Primack 1996). Similarly, 
substructure in galaxy halos can be probed by gravitational lensing even if the satellites contain too few 
baryons to be observable. A lump approximating a singular isothermal sphere (SIS) of one-dimensional 
velocity dispersion a = 10 km s _1 (corresponding to a mass of order 10 8 M Q ) produces a deflection angle 
of about 3 milli-arcsecond if it acts as a lone lens. If the lump resides in a larger halo this deflection will 
effectively be magnified. This can cause distortions of images on ten milli-arcsecond scales and appreciable 
changes to the total flux. 

It has been pointed out by Mao & Schneider (1998) that substructure in the lens galaxy can explain 
the discrepancy between the observed and ('simple') model-predicted flux ratios in the quadruply-imaged 
QSO B1422+231. These authors considered two illustrative pictures of substructure, point-like "globular 
clusters" with masses of 10 6 M Q and a plane-wave smooth fluctuation. In this paper we study how to 
use gravitationally lensed QSOs as a probe of the existence of dark matter subclumps within the lensing 
potential, the small-scale structure that is indeed predicted by numerical CDM simulations. The format is 
as follows. In § ^ we discuss the survival of subclumps in dark matter halos. A model for halo substructure 
is presented in § |^ where it is argued that many subclumps will survive to within a few kpc of a galactic 
halo's center. The formalism for describing lensings by subclumps is introduced in § 0, and the numerical 
simulations used to study the magnification and image distribution are described in § |5j. Specific examples 
of lensing systems and subclump distributions are investigated in § ||. Finally, we discuss the implications 
of this work in § 



2. Survival of substructure in galaxy halos 

In hierarchical cosmogonies - of which CDM is the best studied example - all the mass of the universe 
at early times is contained in small-mass halos that later merge into larger and larger systems. The problem 
of how "halos orbiting within halos" evolve with time is rather complex (Johnston, Spergel, & Hcrnquist 
1995; Moore, Katz, & Lake 1996; Klypin et al. 1999a). Survival of satellites in galaxy clusters has been 
investigated by a number of authors (see Tormen, Diaferio, & Syer 1998; Ghigna et al. 1998) where it has been 
found that somewhere between 10% and 20% of the cluster mass remains in bound subclumps. Simulating 
substructure in galaxy halos is much more computationally demanding. Very high resolution cosmological N- 
body simulations (Klypin et al. 1999b; Moore et al. 1999a) and semi-analytic modeling (Bullock, Kravtsov, 
& Weinberg 2000) have shown that a large number of subclumps do survive in CDM theories, and that 
galaxy halos resemble scaled down versions of galaxy clusters. According to these simulations, the mass 
function of subclumps is approximately dn/dm oc m~ 2 down to a typical resolved mass of ^, 10 7 M . 

To understand the destruction process of dark matter clumps and estimate the range of substructure 
masses and radii that are to be expected, we can use some simple analytic arguments. Dynamical friction 
will cause only the most massive satellites, m £ 10 10 M Q , to sink into the center of galaxy size halos. For 
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the range of (smaller) masses we are interested in, the most important factor in the evolution of substructure 
is tidal stripping within the parent halo that they inhabit. The tidal radius r t of a clump of bound mass m 
can be estimated as 



r t ~R 



3M(R) 



1/3 

(1) 



where R is measured from the center of the (host) halo, and M(R) is the halo mass within that radius 
(Binney & Tremaine 1987). This constraint on the subclump size is plotted in Figure |l| for a Milky Way-like 
galaxy halo approximating a SIS mass density profile, 

P(r) = ^ (2) 

with circular velocity V c = \/2chaio = 226 km s _1 , and for two different radii R. Any clump of radius 
r > rt will be stripped of mass as it falls towards the center. By losing mass and shrinking in size the clump 
may eventually become quasistable when it reaches approximately the tidal stripping line evaluated at the 
pericenter of its orbit (dynamical friction will eventually drag them in on a long timescale and gravitational 
heating caused by the time varying tidal force will slowly pump energy into them). Four such tracks are 
shown in Figure [l] under the assumption that the clump density profile remains unaffected by mass stripping. 

SIS subclumps will lose a large fraction of their mass before they can reach the solar radius, R ~ 10 kpc. 
A clump with a more realistic NFW profile (Navarro, Frenk, & White 1997), 

P{r) = (r/r s )(l+r/r s )^ (3) 

is more centrally concentrated and can reduce its size without loosing much of its mass. This is only true 
until r t (R) ^ r s , at which point the satellite looses its mass more rapidly. Here r s is the subclump scale 
size, i.e. the radius where the logarithmic slope of the profile is —2, p c = 3H 2 /8irG is the critical density, 
and S c indicates the characteristic density contrast of the clump. In Figure [l] the radii r s for three NFW 
tracks are marked with asterisks. If the core of the NFW profile is well below the tidal stripping boundary 
the subclump can survive in the inner halo without losing much mass. According to the r s -S c relation given 
in Navarro et al. (1997), extrapolated to the mass range relevant for this paper (two examples are shown 
in the figure), SCDM halos are right near this boundary of destruction and survival at R ~ 30 kpc. High- 
resolution N-body simulations show, however, a large scatter about this relation (Bullock et al. 1999) and 
in low density models such as currently popular ACDM the halos are more concentrated and thus more able 
to survive - r s is reduced by a factor of ~ 0.6. Since it takes several orbits for a subclump to be completely 
destroyed even if unstable to tidal disruption, and satellites are being continuously accreted by more massive 
galaxies, it is expected that a significant population of subclumps will exist within the Einstein radius of the 
host halo, 

BDl = 4. f^Y ^ « 10 h-' kpc^li^ ( ^ 1 ) 2 , (4) 
V c / D s cD s 1 156 km s- 1 / V ; 

Here Di, D Sl and D\ s are the angular size distances to the lens, source, and from the lens to the source, 

respectively, and Ho = 100 h km s _1 Mpc -1 is the Hubble constant today. The fraction of mass that remains 

in subclumps is investigated more fully in section where a specific substructure model is adopted. 

Another constraint on subclumps comes from the requirement that they not disrupt the observed pop- 
ulation of globular clusters. This is actually more stringent than the constraint imposed by over heating of 
the Galactic disk, which requires a large density of to > 10 6 M Q satellites (Carr & Lacey 1987; Moore & 
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Fig. 1. — Survival of dark matter satellites in a host halo modeled as a singular isothermal sphere of circular velocity 
V c = 226 km B —1 . Dashed lines: tidal stripping radius as a function of subclump mass at two different distances from the 
host halo center, R = 10 and 50 kpc. Dot-dot-dashed lines: constraints imposed by the survival of globular clusters of mass 
m g = 5 X 10 4 Mq and radius r g = 10 pc for 10 Gyr. Dot-dashed line: collisional timescale between subclumps. assuming 
t BB = 10 Gyr (see text for details). Dotted curve: m(r) for a SIS with <r = 5 km s" 1 . Solid curves: three examples of m(r) 
for NFW halos which differ in the choice of r s and <5 C . The asterisks mark r B . The two curves with larger r s follow the r s -8 c 
relation given in Navarro et al. (1997) for the Qm = 1 SCDM model with h = 0.65 extrapolated to these masses. In a ACDM 
model with Hyv = 0.75, r s is about 0.6 times smaller for the same clump mass, so satellites will penetrate further into the halo. 



Silk 1995). If the subclumps are larger than globular clusters, the disruption timescale t sg can be estimated 
in the same way as Spitzer (1958) estimated the timescale for disrupting open star clusters by molecular 
clouds (Binney & Tremaine 1987), i.e. 



E ^ 0.03a ha .iom g r 2 



where r] is the number density of satellites, (Thaio is the velocity dispersion in the parent halo, and m g and 
r g are the mass and radius of the clusters. This constraint is shown in Figure |l] assuming that the fraction 
of the halo mass in satellites is 5%, t sg = 10 Gyr, m g = 5 x 10 4 M©, and r g — 10 pc. It sets a limit on the 
number of dense clumps, not on their size or mass. 

Similarly, the subclump-subclump collisional disruption timescale can be estimated by replacing in 
equation (^|) the mass and radius of a typical globular cluster with those of dark matter substructure, 

t ss — 0.03o-h a io (Gmrijy 1 . (6) 

This is also plotted in Figure |l] for t ss = 10 Gyr. Again, the subclumps do not really need to survive this 
long. Since cosmological simulations show a continuous replenishment of satellites from infall, their number 
could remain high even when collisional disruption is actually taking place. The survival of the baryonic 
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content of infalling clumps is further complicated by shock stripping. Most of the gas may be removed in 
this way, making it difficult to interpret any observations of visible matter in terms of dark satellites, for 
example in the case of the high velocity halo clouds detected in neutral hydrogen (Blitz et al. 1999). 



3. Distribution of satellites 

We seek to develop a simple model for satellites within host halos that incorporates all the aspects of 
real substructure that are important to lensing. We will use a simple power-law model for the subclumps' 
internal structure 

m(m a ) =m a (r/r a ) . (7) 

The parameter r a is kept fixed while the mass within this distance, m a , changes for different clumps. Different 
stages in the stripping of an NFW halo can then be represented by different /3's -(3 = 2 for r t (R) ^ r and 
f3 ~ 0.01 for r t {R) £ r. For the special case of a SIS subclump (3 = 1 and m a = 2a 2 r a /G. Combining this 
with the tidal radius (|l|) gives the mass as a function of galactic radius 



m(m a , R) 



m a I K 



3M(R) \r, 



3=Z 



(8) 



In our model the number density of subclumps, 77, is assumed to be proportional to the average or total 
mass density - clumped plus smooth components. In this way the clumps flow with the mass. Dynamical 
friction is taken to be unimportant. The distribution of satellites is modeled as a power- law in m a , 

1 ""■(">.) = ^(^)^, W 



p dm a 4mor a c 2 \m 

where p is the total local density and m™ ax is the maximum normalized mass. The above is an i?-independent 
quantity. The exponent is chosen so that the true mass function is drj /dm oc m~ n where the normalization 
will depend on R and the structure of the host halo. CDM simulations favor n ~ 2. 

The minimum normalized mass will be set by the parameter f m = m™ ln /w™ ax : the requirement that 
the mass in clumps cannot exceed the total halo mass sets a lower bound on f m if the mass function is 
steep enough. There is a mass below which satellites do not contribute significantly to the lensing and this 
cutoff is in practice more restrictive than the former bound. For this reason the lensing properties become 
insensitive to f m when it is below some level. 

With equations (j|) and (|J) we can easily calculate the mass fraction in subclumps as a function of 
galactic radius for fixed j3. Integrating this along the line of sight gives the fraction of surface density in 
substructure, fs(R). This has been done in figure]^ for a NFW host halo. It can be seen here that for small 
(3, which approximates the case of NFW subclumps with r s £ r t (R), fs(R) remains relatively independent 
of R. By the definition of a subclump r s < R s so we expect that substructure will survive relatively intact in 
the region R R s . As R further decreases (R) will slowly decrees until r t — r s at which point it will drop 
more steeply. It seems realistic that at R ~ 1 kpc, the scale on which multiple QSO images are formed, /s(-R) 
does not drop below 20% of the what it is at R = R s (~ 20 kpc for the Milky Way). Numerical simulations 
find that 10 — 20% of the mass within the virial radius of a 10 12 M halo is contained in substructures of 
mass > 10 7 M Q (Moore et al. 1999a). The lensing properties will be sensitive to structures that are several 
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Fig. 2. — The fraction of surface density in substructure as a function of radius from the center of a NFW halo. The scale size 
of the host halo is r 3 = Rs (see equation ^). This is independent of <5 C . In all cases the normalization is chosen so fs{Rs) = 1. 
SIS subclumps are labeled ,8 = 1. The core of NFW profiles is j3 = 2 and the core of a Moore profile (Moore et al. 1999b) is 
P = 1.5. For comparison R s is expected to be ~ 20 kpc for the Milky Way halo. 



orders of magnitude smaller - below the resolution of any existing simulation. This simple model strongly 
supports our conjecture that in the CDM model /s(-R) is greater than a few percent at the point where 
multiple images are formed and could potentially be much larger. An accurate determination of fs(R) for 
the relevant mass range will require CDM simulations with mass resolution 4—5 orders of magnitude better 
then those available today. 

For the lensing simulations discussed later the subclumps will be modeled as SISs and the distribution 
will be renormalized to the desired /j> The distributions of magnifications and image distortions will be 
only weakly dependent on the SIS assumption because the angular sizes of the sources we consider are not 
small compared to the sizes of the subclumps, as will be seen in section 6 . 2 . 2j . In this case the magnification 
averaged over the image will not depend on the details of the subclumps' internal structure. The lensing is 
more dependent on /s which makes it a more useful free parameter. For SIS subclumps m a can be replaced 
by a in the above giving 



r t (a,R) = a a 



2R 3 



3GM(R) 



m(cr, R) = <t 3 



8R 3 



3G 3 M{R) 



(10) 



1 dr\ 1 

p da m (T max 



2-3n 



(11) 



4. Lensing by substructure 



A probe of purely gravitational mass is needed to confirm or disprove the existence of substructure in 
galaxy halos. We are interested in the situation where there is one galaxy size gravitational lens creating 
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multiplc images of a background source. Within this lens there is a large number of randomly positioned 
subclumps that affect each of the images. We will ignore the effects of any additional clumps outside of the 
main lens. The large number of subclumps and the small, but not vanishingly small, size of the source will 
require simulations of a small region of the lens plane rather than the entire lens at once. The contribution 
from the parts of the lens outside the simulation region must therefore be represented in some way. The 
theoretical background to our method is presented in this section. 

If no substructure existed in the lens it would have the surface density K smoot h(^') expressed in units 
of the critical density, E c = c 2 D s (4:7rGDiD\ s )~ 1 . We will subdivide this density into two components within 
the region of the lens plane being simulated; the substructure component, K su b(a?'); an d the halo component, 

f KsmoothO?') , outside region 

KhaloOr = <^ )_ i . . (12) 

I [1 - .h{x )\ K smoo t h {x ) , mside region 

where x is the center of the region. The simulation region is made large enough that, for the purposes of 
lensing, the mass outside the region can be considered smoothly distributed. It is assumed that f^(x') does 
not change over the simulation region. 

In the thin lens approximation, the light deflection angle caused by the two surface density components 
will add. This follows from linear superposition of the surface potentials. The lensing equation can thus be 
written 

y 1 = x' - a haio (x') - a suh (x'). (13) 

Here the unlensed angular position of the source, 6 S , the position of the image, 9i, and the deflection angle, 
a!{x), have been rescaled 

x = — — y = — — a(x ) = ---a (X x) (14) 

A A A U S 

where A is any arbitrary scale that can be set to a convenient value (usual making <~ 1). In the 

absence of substructure the image in question will appear at xq ■ 

Vo=x — a S mooth(^o)- (15) 
Subtracting equation ( [l5| ) from (|l^), and shifting to the coordinates x = x ' — xq and y = y 1 — yo, yields 

y=X - (3halo(^ + Xq) + <3 smoo th0?o) - a suh {x + Xq). (16) 

We need to express dhaio{x) in terms of a S mooth(a?)- We make the approximation K S mooth(a;') — Khaio(x') = 
f'sK smoo th(x ) inside the simulation region. The constant surface density results in a linear, isotropic deflec- 
tion term fs2K S mooth(x ) x so that the final lensing equation is 

y = X- a smoot h{x + Xq) + d smoot h(x Q ) + fs(x )K(x ) X - d suh (x + Xq). (17) 

The magnification matrix, = dy l /dx^, of equation ( p^ ) is the true observed magnification matrix. This 
is the lensing equation used throughout the rest of this paper. 

One might be concerned that besides the isotropic linear term, /s(xo)k(xo) x, in equation ( |l7j ) there 
might also be anisotropic linear, or shear, term of the same order. This is the shear caused by the surface 
density fsK smoo th(x+xo) inside the simulation region. If the region is small compared to the size of the whole 
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halo this shear turns out to be very small. To see this consider a spherically symmetric mass distribution. 
The contribution to the shear at point R from the mass within the shell R > r > R — 5r is the total 
shear, 7 (r), minus what the shear would be were the shell empty. Since 7(9-) oc r~ 2 outside any spherically 
symmetric mass, and the mass exterior to R does not contribute to j(R), the contribution of this shell is 



Tsholl = j(R) 



(R — Sr) 2 j(R — Sr) 
1 i? 2 7 (i?) 



(18) 



For a singular isothermal sphere j(r) oc r" 1 , so 7s heil/7(-R) = Sr/R. As long as the simulation region is small 
compared to the distance to the center of the host halo the mass exterior to this region will dominate the 
shear in the smooth halo model. If the halo is asymmetric it is only asymmetries on a scale smaller than the 
simulation region that could change this conclusion. In all the cases we will consider, the simulation region 
will be a factor of 100 or more smaller than the distance between the image and the center of the lens, R. 
We can then safely ignore this shear term as is done in equation (|l7|). Another related consideration is that, 
in the limit where the number of subclumps becomes very large at constant surface density, the distortion of 
the image should converge to the smooth lens solution. In our simulations the satellites will be distributed 
uniformly over the region considered so it is incapable of producing a shear in this limit. As a result the 
extra shear must be zero. 

It is instructive at this point to estimate how large of an effect substructure is expected to have. 
Expanding equation (|l7| ) around xq we can estimate the distortion of an image along the two eigenvectors 
of the smooth model's magnification matrix, 

91 ~ l-(l-.fc)«±7- (19) 

It is apparent that the effects of substructure will be magnified or demagnified by the larger scale lens as 
represented by the denominator in equation (|l9|) . Dark matter subclumps will have a significant effect on the 
image only when the deflection is a significant fraction of the angular size of the source 9 S £ (a S ub)A /Dz. 
The question of whether the lensing effect of satellites is detectable reduces then to the question of whether 
deflections of this order are common. Note that the probability of a subclump contributing to the lensing 
is boosted by a factor equal to the background magnification /i b = ([1 — (1 — ./s)k] 2 — 7 2 ) , which can be 
very large near a critical curve where the magnification formally diverges. This, combined with the fact that 
the line of sight already passes through a high density region, result in substructure having a larger effect 
than might otherwise be expected. 



5. Lensing simulations 

To calculate the distribution of magnifications for a large number of random configurations of satellites 
an efficient method of simulating the lensing effect is required. The basic approach we use is to calculate the 
deflection angle, a su b(a?), at every point on a dense grid, and use equation ( |l7| ) to determine if y lies within 
a specified source. Since surface brightness is conserved the magnification can be calculated by multiplying 
the area of a grid cell by the number of grid points found to be in images and dividing by the original size 
of the source. We are considering images that are smaller than the angular size of the host halo, so the grid 
will cover a small region over which the average density is approximated as constant. Multiple sources can 
be considered simultaneously as long as it is certain that none of their images are outside of the gridded 
region. 
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Direct summation of the contributions to a su b(x) from each subclump is ponderously slow (of order 
-^ciump^Vgrid)- To speed the calculation up the following approach is used. We scale all lengths to a convenient 
and fixed A Q . After a random realization of subclump masses and positions is created, a su b(%) is calculated 
in two steps. The mass distribution in clumps can be related to the deflection angle through the lcnsing 
potential, which can be broken into two parts %j)(x) = 5ip(x) + i/)°(x). The background potential, ip°(x), is 
the solution to V 2 ^°(x) —2k where k is the average surface density of substructure within the simulation 
region. This is easily solved. The perturbed potential and the deflection angle are then given by 

V 2 5ip(x) = 2 [k(x) - k] , a suh (x) = V8tp(x) + kx. (20) 

An estimate of a su b(^) can be quickly found by moving the mass of each clump to its closest grid point and 



solving (20) by fast Fourier transform (FFT). 



The contribution to a su b(%) from nearby subclumps is not well represented in this FFT solution because 
of the interpolation to a grid and the internal structure of a subclump. To refine the calculation a square 
patch around each clump is considered. The contribution of this satellite to the FFT estimate is subtracted 
in this patch, and then a more accurate estimate is recalculated. Here we use the fact that the deflection 
angle for any spherically symmetric lens is particularly simple (Schneider, Ehlers, & Falco 1992) 

a{ x ) = x m{x) =2 f dx'x'n{x'). (21) 



•*> Jo 

The size of the patch around each satellite is adaptively scaled to be proportional to the size of the clump. 
This approach results in a ~ N gI id log A gr id + A c i ump A patc h scaling where A patc h is the average number of 
grid points in the patches. Generally there are more small subclumps than large ones, and this technique is 
much faster than the direct sum over satellites. 

In this way the FFT is used to calculate the long range "forces" , and direct summation is used for short 
range as in a P3M N-body code. The FFT results in the effective mass distribution being periodic on the 
scale of the box; in practice, however, the mass distribution is always smooth enough on these scales and 
the artificial periodicity does not affect the lensing. 

Note that if k = /£Ksmooth(£o), as it is on average, the second term in the expression for a S ub(af) in ( ^p| ) 



will cancel the fourth term in equation (17). This will generally be the case when there are large number of 
subclumps within the simulation region making statistical fluctuations small. This ensures that the results 
will converge to the smooth lens case when the subclumps get very dense. In the opposite limit where there 
are no subclumps in the simulation region a su b(x) — and equation ( |l7| ) will give the correct background 
magnification. 



6. Examples 

A comparison of our simulations with the data requires a specific lensing system to be modeled in detail. 
An overall smooth lensing potential needs to be constructed before the effects of substructure on each image 
can be investigated. We will leave a detailed analysis of specific lensing systems for a later paper, and concern 
ourselves here with representative examples that will illustrate what can be expected. In section 3.1 some 



general expectations for the SIS model that is used in the simulations will be discussed. Then in section 
the results of numerical simulations are presented and interpreted. 



6.2 
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6.1. SIS model 

For simplicity, the host halo will be modeled as a SIS. While to model most observed lens systems some 
ellipticity or external shear must be added to the SIS model, the effects of substructure will not qualitatively 
change in an asymmetric model. 

Two images will then appear at x = y±l (y < 1), where the scale used is A (c7haio) = 47rcrf; alo c~ 2 Z?iZ?i s /.D s . 
The unperturbed shear and convergence are 

2r 2\x\ 

From equation (|l0|), the subclump - also SISs - have sizes and masses 

r(a, R) = m(a, R) = ™f . (23) 

V O (T halo V oGcr halo 

The scaled deflection angle caused by a truncated SIS subclump can be calculated using equation (|2l|), 
„, s 2 / a \ 2 ( a - Va 2 - 1 + tan -1 V a 2 - 1 a > 1 

a \ x ) = - - — M _ _ ^ -i ( 24 ) 



a , a < 1 

where a = r(a, R)/xX (tj max ), and x is the distance defined in equation ( fTi] ) using A(<7 max ). 

When the impact parameter y is small the images will form at close to the same R on opposite sides 
of the lens. In section || it was shown that if the subclumps are compact enough compared to the host halo 
their mass distribution will not be a strong function of R. Following this argument we will take fy, to be 
the same for both images in the simulations. We will also replace the mass and radius distributions at a 
projected distance x with their values at the radius R = x\ {a\ iai \ ) . Most of the mass to located near this 
radius and a more detailed model would not be warranted at this point. 

Let us consider what level of substructure will be important for lensing in this model. From § |i] we know 
subclumps will have a significant effect on the image when the deflection caused by an individual satellite 
is a significant fraction of the angular size of the source. We will use the requirement 8 S ^ 10 (a su b) A / D\. 
A characteristic deflection angle is (a su b)A — A (er). We can then put a loose lower limit on the mass of 
subclumps that can have a significant effect on the lensed image of a source of physical size l s , 
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where the tidal radius ( p3| ) was used. As we shall see this estimate is fairly accurate. 

If subclumps with m > m c are rare because the mass fraction in this mass range is small, then it will 
be unlikely for images to be affected. We can introduce a figure of merit that gauges the effectiveness of 
substructure to influence images by considering the average number of satellites that intersect an image. An 
idealized circular image will have a radius of l s Diy/\ij\/ D s on the lens plane. If we add this radius to the 
radius of the subclump we can find the area in which a subclump will influence the image. The average 
number of subclumps in this region is, 



F mi u = 7r / dm' 
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Fig. 3. — Some examples of the substructure lensing parameter -F su b as a function of source size. The solid line 
is for (/e; Mi R, mmai) = (0.05, 9.3, 1 kpc, 10 10 Mq) where /s in this case is the fraction of mass in subclumps with 
m > m c . Parameterized in the same way - dotted (0.1,7.3,1 kpc, 10 9 Mq), dashed (0.03,5.0,1 kpc, 10 9 M ), dot-dashed 
(0.03, 3.0, 5 kpc, 10 9 Mq), dot-dot-dashed (0.05, 9.3, 1 kpc, 10 s Mq). In all cases z s = 3, 2; = 1, the subclump mass function is 
oc m~ 2 , and we assumed a flat cosmology with Qjvf = 1- 

Roughly, if -F su b ^ 1 substructure will affect almost every image. How well ~F su b characterizes the lensing 
properties needs to be investigated with numerical simulations, but we can glean some expectations from 
Figure || where F su b is plotted for some example sets of satellite parameters. Clearly, only a small fraction of 
the mass needs be in subclumps, ~ 1%, for them to have a significant influence on images. Generally, when 
the source is small substructure is important because m c is small and there are many available lenses. As 
l s increases -FLub decreases until the source gets large enough that it becomes probable for one of the larger 
satellites to lie close to the image. F BU b would not increase again for large l s as it does in Figure | if the mass 
function were more bottom-heavy than mT 2 . When the grow subclumps larger while keeping their masses 
fixed the minimum of -F su b is smaller and moves to larger l s \ the parameter R is used here to normalize the 
size distribution through equation (p3|). -FLub does not increase indefinitely with l s because at some point 
w max ~ m c (/ s ), and satellites cease to have much of an effect. 

The source has been taken to be circular so far, but if the source is highly elongated like a radio jet 
the likelihood of there being a subclump nearby will be larger. At the same time its small width makes it 
susceptible to distortion by comparatively small subclumps. The superior resolution of radio interferometers 
could make these good sources for observing substructure lensing if the surface brightness of the jets is high 
enough. 

Microlensing of multiple image QSOs has already been identified through its time dependence. Stellar 
mass compact objects have lensing timescales on the order of months. In contrast lensing by substructure 
will be nearly time independent. The characteristic timescale in this case is t = A (er)/uj_. This results in 

r , / a \2 / kmr^ /H DiD ls \ 
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Fig. 4. — The cumulative magnification distributions for the two images of a SIS lens with source position y = 0.12 measured 
in lensing scale lengths. The distribution of image 1 - the primary image or outer image - is plotted on the left. In the center 
is the distribution for image 2. On the right the cumulative distribution of the difference in image brightnesses measured in 
magnitudes is shown. Negative values of Am correspond to cases where image 2 is brighter than image 1. The vertical dotted 
lines are the values expected for a smooth lens and an infinitely small source. The different curves are for different forms of 
substructure and source sizes. The solid curve assumes that 5% of the mass is in subclumps with 10 4 Mq < m < 10 8 Mq, and 
a 10 pc source. Dashed line: same except with 10% of the mass in substructure. Dash-dotted curve: assumes 5% of the mass 
in subclumps with 10 3 Mq < m < 10 7 Mq and a 1 pc source. Dotted curve: same as the dash-dotted line only with 10% of 
the mass in substructure. 



For subclumps with cr su b ~ 10 km s 1 and v± ~ 200 km s 1 this timescale is very large. It is therefore 
justified to treat the subclump as stationary. 



6.2. Numerical results 



In the following, different aspects of our simulation results will be discussed which correspond to different 
methods for probing substructure. The simulations are carried out for a compact (< 1 kpc) source at z = 3 
in an Einstein-de Sitter cosmology. The cosmology has little effect on the results. The subclump mass 
function is cx rn~ 2 in accordance with expectation from CDM models. A 256 2 grid is usually used giving 
an estimated numerical error in the magnification of <5/i ~ ^4. g rid(-^ 2 A sourco ) _1 = 10~ 5 A gr id/A sourcc where 
the A's are the area of the grid and source. The grid area changes according to the subclump sizes and 
the background magnification, but generally an accuracy of Sfi ~ 0.01 — 0.05 is attained. Simulations with 
higher resolution have been run to verify this accuracy estimate. 



6.2.1. Magnifications 



As will be seen in section 6.2.2 the image positions are not significantly changed by small scale sub- 
structure. This allows for the possibility of constructing a smooth lens model from the image positions and 
perhaps the images of the QSO's host galaxy. The magnification ratios predicted by this model, or family of 
models, can then be compared with the observations. This will generally require four image systems because 
two image systems do not contain enough information to constrain the smooth lens model. For simplicity 
we will investigate a simple SIS lens at z = 1 with a velocity dispersion of Chaio = 240 km s ; this is a high 



- 13 - 




6 8 10 12 14 2 4 6 8 10 12 -1.0-0.5 0.0 0.5 1.0 1.5 2.0 

\fj.\ Am 



Fig. 6. — The cumulative magnification distributions for the two images of a 10 pc source with source position y = 0.12. 
The three different curves correspond to different subclump mass ranges. Solid line: 10 4 Mq < m < 10 5 Mq. Dashed line: 
10 6 M < m < 10 7 Mq. Dotted line: 10 7 M < m < 10 8 M a . The mass fraction in subclumps, / s , is held fixed at 2%. The 
subclumps are progressively more influential on the lensing as they get more massive despite their smaller number density. The 
vertical sections of the dotted curves are located at the background magnifications - the value of the magnification when there 
are no subclumps influencing the image. 

velocity dispersion for galaxies in general, but the population of gravitational lenses tend to be highly biased 
toward the largest of elliptical galaxies. In any case, if the source position is fixed in units of A (tTh a io) then 
Chaio only affects the subclump size distribution - a smaller host halo giving smaller subclumps at the image 
positions. 

The magnification distribution for each image is calculated by repeating the calculation with different 
realizations of the substructure. Figures ^ and ^ show the cumulative distributions for two choices of the 
projected source position relative to the center of the host halo. The magnification ratios of the images 
- the only directly observable quantity - are shown on the right of each figure. This demonstrates that 
with only a few percent of the surface density contained in substructure the deviations from the expected 
magnification ratio are significant in most cases. For comparison typical errors quoted for HST observations 
are ^ 0.05 mag. For the models in Figures [| only 14% — 32% of the cases, depending on the model, have 
ratios within 0.2 mag of the expected ratio. For Figures || where the source position is further from the lens 
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Fig. 7. — Same as in Figure pi but the source is 1 pc in size and the mass ranges are: Solid line: 10 2 Mq < m < 10 3 Mq. 
Dashed line: 10 4 Mq < m < lCr Mq. Dotted line: 10 5 M < m < 10 6 Mq. 



center this is 26% — 44% of the cases. For the models in Figure [| a significant fraction of the cases have the 
image that is expected to be the brightest instead the dimmer of the two - a sure sign that the lens is more 
complicated than the lens model reflects. 

To gauge the impact of substructure on lensing data consider n lens systems. The likelihood that at 
least one of the magnification ratios is discrepant by a certain amount, m, is 1 — P1P2 ■ ■ - Pn where pi is the 
probability that the ith image pair has a magnification ratio within m of the expected ratio. With only a 
few systems a strong limit can be put on the amount of substructure as long as the expected magnification 
ratios can be predicted to a certain degree. For example, if we were to observe 5 independent QSO image 
pairs like the ones depicted in Figures ^ or |j| and all of them had ratios within 0.2 mag of the expected 
values then all the models displayed could be ruled out at better than 98% confidence. 

We expect that if l\ / m 2 is held fixed the magnification distribution will be relatively independent of the 



source size and the subclump mass range; this is because A (cr) cx a oc m ' (see § 6.1). A comparison of 
figures |^ and || is generally in agreement with this scaling relation. Without precise knowledge of the source 
size the strongest constraint would be put on fa. If the substructure masses are reduced while keeping the 
source size fixed the distribution tends not to extend to as high magnifications because there is an effective 
limit to how much a small structure can magnify the large sources. For a fixed substructure distribution, as 
expected from our discussion of F su b (section |6~l| ) , the images of smaller sources tend to be more strongly 
affected by substructure. Smaller sources are sensitive to a wider range of subclump masses. Figures |^ and |7] 
show how the lensing depends on subclump mass. Subclumps of mass 10 4 M Q do not have a significant 
influence on 10 pc sources and m ^ 10 2 M subclumps do not affect 1 pc sources at the position and redshift 
considered. This is in general agreement with the estimate of m c made before, equation (|2^). This is also a 
useful test of the simulations because as the number density of subclumps goes up the magnification should 
converge to the smooth lens solution. 

For the high substructure mass ranges in Figure ^ and especially Figure ^ there is a bias toward smaller 
than expected magnification ratios. This is a result of substructure having a different effect on different 
kinds of images. To understand this consider a small change in the convergence, k, (surface density) from 
the value predicted by the smooth model k d . To first order the change in the magnification is 

— ^ ~ 2(1 - k )(i 5k . (28) 
Mo 
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Clumping increases the probability that an image will land on an underdence region where 6 k < 0. The 
primary image - where fi a > 1 and k < 1 - will have a smaller magnification in this case. For the secondary 
image, which is reflected with respect to the primary, the situation is more complicated. If y < 0.5, like it is in 
our examples, then ji < and k < 1, which means the absolute value of the magnification will go up when 
K goes down. The opposite is true when y > 0.5 since in this case fi < and n a > 1 for image two. This 
biasing when y < 0.5 does not seem to have a large effect when the mass distribution is steep and extends 
down to m c as in Figures ^| and ^. However, if the substructure consists of only high mass objects - relative 
to the source size - the bias will be significant. This is most noticeable for the 10 7 M© < m < 10 8 Mq 
model in Figure || (dotted curves) where in a large fraction of the cases the images have the background 
magnification, fi b (the vertical sections of the curves). 

When confronting observations the smooth lens model and the substructure model need to be constrained 
simultaneously in four image lens systems. This task is complicated by degeneracies in smooth lens models 
which produce equally good fits to the image positions, but different magnification ratios. Constraining the 
level of substructure in these systems requires significant further development beyond what is presented here 
and will be the subject of a future paper by one of the authors (Metcalf & Zhao 2001). 



6.2.2. Image distortions 

Resolution in the visible is not high enough to directly observe distortions to images caused by substruc- 
ture in the mass range that we are considering here. However radio observations with sub milli-arcsecond 
resolution could be used to find distortions in one image that are not mirrored in all the other images if 
the source is of the appropriate size. To investigate this possibility and to gain more understanding of the 
lensing behavior in general we consider the images produced in the same simulations used in the previous 
subsection. 

Figure ^ shows representations of the images with random realizations of substructure. In most cases 
the change in the morphology of the image is small. The distortions are typically on milli-arcsecond scales. 
It can also be seen that the centroid of the image is not significantly shifted. It would take substructures 
with masses of <L 10 8 M that are very well aligned with the image to change the image position by a few 
tens of milli-arcsecond. Such a chance alignment will be rare even in the CDM model. 

These figures also show that with the sources and substructures we are considering the internal structure 
of the subclumps should not make a great deal of difference to the lensing. The images feel a deflection 
potential that is smoothed on the scale of their own size. Since the images are about the same size as the 
subclumps in these examples how cuspy the subclump cores are will not make a great deal of difference. 

Another interesting possibility is a lensed radio jet. Jets typically have opening angles of a few degrees 
and extend over Mpc (see Kembhavi & Narlika 1999 and Peterson 1997). This source would then cover more 
area on the lens plane making it more likely to come close to a large subclump; its thinness would make 
distortions more pronounced. Figure [)] shows some idealized representations of what a lensed radio jet would 
look like. Since the jets are not perfectly conical in reality the distortion caused by lensing would need to be 
deduced by comparing the two images. Some of the distortions are on ~ 0.01 arcsec scales. There is at least 
one multiply-imaged radio jet known, Q0957 + 561. The two images have been mapped on milli-arcsecond 
scales (Gorenstein et al. 1988) and there is not any obvious evidence for substructure in this case although 
a rather large subclump would be require to produce a noticeable effect (~ 10 7 M Q )(Barkana et al. 1999). 
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6.2.3. Distribution of magnification ratios 

Another approach to take in the search for substructure is to look at all multiply-imaged QSOs and con- 
sider the overall distribution of magnification ratios without modeling individual cases, only the population 
of halos in general. This method has several potential advantages besides requiring less work. Two image 
QSOs can be used in this case, which increases the size of the sample and generally the physical distance the 
images are from the center of the lens. The latter property means both that substructure will be more likely 
to survive and that the final constraint will be less dependent on the smooth lens model used. To investigate 
this situation we fix the source at z = 3, but choose a random lens redshift and impact parameter according 
to the assumption of constant comoving density. The physical size distribution of subclumps is kept fixed 
and is the same as in the previous two subsections. In this case the lensing distribution is independent of 
the distribution of lens velocity dispersions. 

The cumulative distribution of magnification ratios is shown in Figure [l0| along with the theoretical 
prediction for smooth lenses. The most significant effect is an over abundance of cases with magnification 
ratios close to one and cases where the image that is expected to be brighter is in fact the dimmer image. 
This is mostly the result of simple spillover caused by extra scatter added to a steep distribution - the biasing 



discussed in section |6.2.l| does not play a big part. Identifying the image that is expected to be brighter 
requires some knowledge of the individual lenses. In the simple case of an SIS lens it would just require 
determining which image is closer to the lens; but without this information - which is often missing in radio 
surveys - we can still consider the ratio of the dimmer image to the brighter one. This distribution is the 



dashed one in Figure 10 where it can be seen that there is an overabundance of cases with ratios close to 
one. A proper comparison with data would require accurate modeling of the selection effects inherent in any 
lens survey as well as using a more realistic smooth lens model. At this time the publicly available data is 
very limited and does not put an interesting constraint on substructure. This may change in the near future 
with the completion of several large, systematic QSO surveys. 



7. Discussion 

We have shown in this Paper that, if dark matter substructure exists in the halos of galaxies and can 
survive to constitute a percent or more of the surface density at small impact parameters (several kpc) - a 
possibility that is in general agreement with CDM hierarchical cosmologies, it will cause significant changes 
to the magnification ratios of multiply-imaged QSOs. Furthermore, if the magnification ratios of just a few 
QSO image pairs can be predicted from the image positions and other information to an accuracy of a few 
tenths of a magnitude the existence of such substructure can be confirmed or invalidated. This would allow 
us to probe primordial structure formation on the scales of 10 3 — 10 8 M Q . To confront observations more 
detailed modeling of individual QSO lens systems will be needed. This is the topic of another paper (Metcalf 
& Zhao 2001, in preparation). 

Besides lens modeling there are several other complicating factors involved in measuring the magnifi- 
cation ratios. A correction for differential extinction must be included in the visible bands. Falco et al. 
(1999) find a median differential extinction between lensed images of AE(U — V) = 0.04 — 0.06 mag so this 
is not likely to be a large impediment. Another issue is intrinsic variation in the source. If the source is 
variable and the time delay is not well known it is difficult to interpret the flux ratios. Radio observations 
avoid the complication of differential extinction, but radio galaxies are often variable. Also the lens galaxy is 
usually not detected in the radio, making the constraints on the smooth lens model weaker. A combination 
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of observations at different wavelengths will be required. 

Microlensing by stars in the lens galaxy is also a possible contaminating factor if the source is too 
small. The continuum emission region of a QSO is believed to be less then ~ 100 au, and is observed to be 
microlensed in some cases. To avoid this a larger source must be used. The broad line emission region is 
believed to be about a pc in size, which would make it sensitive the subclump masses £ 10 3 M . There are 
also features in the radio emission that are both smaller and larger than this. They would generally need to 
be resolved in order to confirm their scale. QSO narrow line regions probably have sizes of ^ 100 pc. The 
molecular line emission region of AGN is believed to be 1 — 100 kpc in size, but it is not certain if the CO 
emission from QSOs, for example, comes from a localized region or from the host galaxy as a whole (see 
Kembhavi & Narlika 1999 for a review of QSO properties). Conclusions drawn from lensing on the existence 
and mass distribution of substructure will thus depend on the kind of observations used and our knowledge 
of the sources' internal structure. 

Finally, the images could also be lensed by more familiar galactic substructures such as spiral arms (Mao 
& Schneider 1998). These are expected to interfere with the constraints on dark matter structure at a small 
level. The Milky Way disk has a surface density of ~ 50 Mq/ pc 2 and spiral arms are believed to have an 
amplitude of about 20% in mass. For a lens at z — 1 and a source at z = 3 in an Einstein-de Sitter universe 
this corresponds to ft = 0.003/i _1 . We can use equation ( p8| ) to calculate the change this would cause in the 
magnification, which is ~ 0.05 mag in the case of the most magnified image discussed in section ^. The dark 
matter subclumps considered here have more influence on the magnification because of their comparatively 
high mass density. Globular clusters are unlikely to cause significant changes in the magnifications because 
they do not contain a significant fraction of the halo mass (~ 10~ 4 of the Milky Way halo) as pointed out 
by Mao & Schneider (1998). 
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Fig. 8. — Images of a 10 pc source at z = 3 seen through a lens at z = 1. On the left are realizations of the outer, primary 
image and on the right simulations of the inner image. The source position is y = 0.12. The subclumps constitute 5% of the 
density and 10 4 M© < m < 10 s Mq with a m~ 2 mass function. For clarity only a small part of the field used in the simulations 
is shown and the subclumps with r < 0.3\ o (m max ) are not drawn. This corresponds to masses below 2.0 X 10 4 Mq. The radii 
of the circles are the tidal radii of the subclumps. The thick dashed curves are the shape the image would be if the lens were 
smooth and the thick circles arc what the source would look like if there were no lensing. 
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Fig. 9. — Images of an idealized radio jet at z = 3 seen through a lens at z = 1. The source of the jet is at y = 0.12. 
The subclumps constitute 5% of the density and 10 4 Mq < m < 10 8 Mq with a mT 2 mass function. Only subclumps with 
r > 0.35A o (m 

max 

) (m > 3.2 x 10 4 M Q ) are shown. The unlensed jet is shown in dotted lines. The opening angle is 6° and it 
is oriented at 10° from the lens' radial direction that is the vertical axis here. The dashed lines show the image were the lens 
smooth. 
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Fig. 10. — The distribution of magnification ratios for a 10 pc source at z = 3. The smooth curve is for lenses without 
substructure. The solid histogram includes the effect of subclumps. The dashed histogram is the same except that Q is the 
lesser of /ii//i2 and /U2//U1. The subclumps are a fixed 5% of the lens surface density, and have 10 4 Mq < m < 10* Mq. All 
three distributions have the same integral for Q > 2/7. This cutoff is necessary because the Monte Carlo calculation cannot be 
efficiently extended down to Q = 0. However, small Q pairs are strongly selected against in observations and the distribution 
is not strongly affected by substructure in this region. There are 10,000 realizations in the simulated distribution. 



